#include <iostream>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <queue>

#include "SparseMatrix.h"

#define REP(i, n) for (int i = 0; i < (n); ++i)
#define FOR(i, a, b) for (int i = (a); i <= (b); ++i)
#define FOREACH(i, coll) for (typeof(coll.begin()) i = coll.begin(); i != coll.end(); ++i)

template <typename GeneratorType, typename ScalarType>
class ReducibleFreeChainComplex {
	typedef int IndexType;

	vector<vector<GeneratorType> > generators_;
	vector<map<GeneratorType, IndexType> > indices_;

	vector<SparseMatrix<IndexType, ScalarType> > bd_;

	int embeddingDimension() const {
		return (int)bd_.size() - 2;
	}
	
	bool find_reduction_pair(int i, IndexType* a, IndexType* b) {
		return bd_[i].find_min(a, b);
	}

	void reduce(int i, const IndexType& a, const IndexType& b) {
//		cout << "reduce(" << a << ", " << b << ")" << endl;

		bd_[i+1].delete_row(b);

		ScalarType coeff = bd_[i].get(a, b);

		map<IndexType, ScalarType> row = bd_[i].get_row(a);
		FOREACH(e_, row) {
			const IndexType& e = e_->first;
			if (e != b) {
				const ScalarType& coeff_e = e_->second;
				bd_[i].sub_col(e, b, coeff * coeff_e); 
			}		
		}
		indices_[i].erase(generators_[i][b]);
		indices_[i-1].erase(generators_[i-1][a]);
		bd_[i].delete_col(b);
		bd_[i-1].delete_col(a);
	}

	public:
	
	ReducibleFreeChainComplex(int embDim, const vector<vector<GeneratorType> >& generators) : generators_(generators) {
		generators_.resize(embDim + 2);
		indices_.resize(embDim + 2);
		REP(i, generators_.size())
			REP(g_i, generators_[i].size())
				indices_[i][generators_[i][g_i]] = g_i;
		bd_.push_back(SparseMatrix<IndexType, ScalarType>(0, generators_[0].size()));
		FOR(i, 1, embDim + 1)
			bd_.push_back(SparseMatrix<IndexType, ScalarType>(generators_[i-1].size(), generators_[i].size()));
	}

	void reduce() {
//		cout << *this << endl;
		for (int i = embeddingDimension(); i >= 1; --i) {
			IndexType a, b;
			while (find_reduction_pair(i, &a, &b)) {
				reduce(i, a, b);
				cout << *this << endl;
			}
		}
	}

	void add_mapping(int i, const GeneratorType& a, const GeneratorType& b, const ScalarType& c) {
		bd_[i].set(indices_[i-1][a], indices_[i][b], c);
	}

	int size(int i) const {
		return indices_[i].size();
	}

	friend ostream& operator<<(ostream& out, ReducibleFreeChainComplex& E) {
		cout << "Chain Complex: " << endl;
		for (int i = 0; i < E.indices_.size(); ++i) if (E.indices_[i].size()) {
			out << i << ": ";
			FOREACH(g_, E.indices_[i]) {
				out << g_->first;
				if (g_ != (--E.indices_[i].end()))
					out << ", ";
			}
			out << endl; 
		}
		for (int i = 0; i < E.bd_.size(); ++i) {
			out << i << ": ";
			REP(g1, E.bd_[i].cols_.size()) if (E.bd_[i].cols_[g1].size()) {
				out << E.generators_[i][g1] << " = ";
				FOREACH(g2, E.bd_[i].cols_[g1]) {
					if (g2->second < 0)
						out << "-";
					else if (g2 != E.bd_[i].cols_[g1].begin())
						out << "+";
					if (abs(g2->second) != 1)
						out << abs(g2->second);
					out << E.generators_[i-1][g2->first];
					typename map<IndexType, ScalarType>::const_iterator iter = g2;
					iter++;
				}
				out << "; ";
			}
			out << endl; 
		}	
		return out;
	} 
};

